In [6]:
from matplotlib import pyplot as plt
plt.imshow(mask)
Out[6]:
<matplotlib.image.AxesImage at 0x2eb7e7bdd00>
There are two quality layers associated with L2A Reflectance Product contained within the L2A Mask file. The mask variable contains 6 binary flag bands (1-5, 8) which should be excluded from analysis, and two data bands (6,7).
The second band_mask variable indicates whether or not any given wavelength of any given pixel is interpolated. Interpolation occurs either due to a focal plane array bad pixel, or from saturation. This data is provided as a packed unsigned integer array with 36 elements.
Note: This How-to guide currently only explains and utilizes the
maskdata, not theband_mask.
Requirements:
emit_tutorials environment as the kernel for this notebook.../data/ folder.Learning Objectives
# Import Packages
import os
import netCDF4 as nc
from osgeo import gdal
import numpy as np
import xarray as xr
import hvplot.xarray
import holoviews as hv
import sys
sys.path.append('../modules/')
import emit_tools
In this guide we will use an EMIT L2A reflectance file as well as an EMIT L2A Mask file. Set the filepaths for each.
fp = '../data/EMIT_L2A_RFL_001_20220903T163129_2224611_012.nc'
fp_mask = '../data/EMIT_L2A_MASK_001_20220903T163129_2224611_012.nc'
The most efficient way to utilize the mask is to apply it before orthorectification because the orthorectified datasets take up more space. To apply a mask using the L2A Mask file, we want to open it, specify which bands to use in construction of a mask, and then apply the mask.
To do this, first take a look at what each band will mask by reading in the sensor_band_parameters group from the mask file as an xarray.dataset then converting to a dataframe.
Note: In the user guide, the bands are indexed as 1-8 not 0-7 as used here.
mask_parameters_ds = xr.open_dataset(fp_mask,engine = 'h5netcdf', group='sensor_band_parameters')
mask_key = mask_parameters_ds['mask_bands'].to_dataframe()
mask_key
| mask_bands | |
|---|---|
| bands | |
| 0 | Cloud flag |
| 1 | Cirrus flag |
| 2 | Water flag |
| 3 | Spacecraft Flag |
| 4 | Dilated Cloud Flag |
| 5 | AOD550 |
| 6 | H2O (g cm-2) |
| 7 | Aggregate Flag |
The above dataframe shows exactly what each band contained within the file represents/will mask. For this example we will use flags 0,1,3, and 4 to remove any potential clouds and any artefacts caused by the space station. This can be done using the build_mask function from the emit_tools module. This function combines the requested flags into a single mask and returns it as an array.
Select the bands to use.
flags = [0,1,3,4]
flags
[0, 1, 3, 4]
Now create the mask using the build_mask function.
mask = emit_tools.build_mask(fp_mask,flags)
Flags used: ['Cloud flag' 'Cirrus flag' 'Spacecraft Flag' 'Dilated Cloud Flag']
To see the regions of the unorthocorrected image that will be masked we can plot the mask array.
from matplotlib import pyplot as plt
plt.imshow(mask)
<matplotlib.image.AxesImage at 0x2eb7e7bdd00>
Now that we have a mask to apply, we can use it as the quality_mask parameter in the emit_xarray function.
Note: If you do not orthocorrect the image, the mask will simply be added as a
quality_maskvariable.
ds = emit_tools.emit_xarray(fp, ortho=True, quality_mask=mask)
ds
<xarray.Dataset>
Dimensions: (latitude: 2009, longitude: 2353, bands: 285)
Coordinates:
* latitude (latitude) float64 -39.31 -39.31 -39.31 ... -40.39 -40.4 -40.4
* longitude (longitude) float64 -62.51 -62.51 -62.51 ... -61.24 -61.24
wavelengths (bands) float32 381.0 388.4 395.8 ... 2.486e+03 2.493e+03
fwhm (bands) float32 8.415 8.415 8.415 8.415 ... 8.806 8.807 8.809
spatial_ref int32 0
Dimensions without coordinates: bands
Data variables:
reflectance (latitude, longitude, bands) float32 nan nan nan ... nan nan
Attributes: (12/38)
ncei_template_version: NCEI_NetCDF_Swath_Template_v2.0
summary: The Earth Surface Mineral Dust Source ...
keywords: Imaging Spectroscopy, minerals, EMIT, ...
Conventions: CF-1.63
sensor: EMIT (Earth Surface Mineral Dust Sourc...
instrument: EMIT
... ...
southernmost_latitude: -40.39610428069674
spatialResolution: 0.000542232520256367
spatial_ref: GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
geotransform: [-6.25120945e+01 5.42232520e-04 -0.00...
day_night_flag: Day
title: EMIT L2A Surface Reflectance 60 m V001This ds dataset is now orthorectified and the specified mask has been applied. Visualize the output using a plot of a red band (650 nm).
red = np.nanargmin(abs(ds['wavelengths'].data-650)) # Find closes band to 650 - Red
ds.isel(bands=red).hvplot.image(cmap='viridis', aspect = 'equal', frame_width=500, rasterize=True)
The modified xarray.Dataset can also be saved as a netCDF4 output that can be reopened using the xarray.open_dataset function.
ds.to_netcdf('../data/example_quality_nc_out.nc')
# Example for Opening
# ds = xr.open_dataset('../data/example_quality_nc_out.nc')
Email: LPDAAC@usgs.gov
Voice: +1-866-573-3222
Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹
Website: https://lpdaac.usgs.gov/
Date last modified: 01-04-2023
¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.